using System;
using MetaViewAnisotropicDiffusion.Matrix;
using MetaViewAnisotropicDiffusion.Utility;
using MatrixDataStructures;

namespace MetaViewAnisotropicDiffusion
{
    namespace PeronaMalik
    {
        public class Anisotropic
        {

            public Anisotropic(Single dt)
            {
                lambda = dt;
            }

            protected static double lambda = 0.25;

# region "For GrayScale Bitmap 2D"

            public static void AnisoDiff(DoubleMatrix image, double K, int iterations)
            {
                for (int i = 0; i < iterations; i++)
                {
                    AnisoDiff(image, K);
                }

            }
            public static void AnisoDiff(DoubleMatrix image, double K)
            {
                ComplexDoubleMatrix grad = new ComplexDoubleMatrix(image.XDim, image.YDim, 0d);
                ComplexDoubleMatrix flux = new ComplexDoubleMatrix(image.XDim, image.YDim, 0d);

                // Calculate gradient
                for (int i = 1; i < image.XDim; i++)
                {
                    for (int j = 0; j < image.YDim; j++)
                    {
                        grad.real.data[i, j] = image.data[i - 1, j] - image.data[i, j];
                    }
                }

                for (int i = 0; i < image.XDim; i++)
                {
                    for (int j = 1; j < image.YDim; j++)
                    {
                        grad.imag.data[i, j] = image.data[i, j - 1] - image.data[i, j];
                    }
                }


                // Calculate flux
                for (int i = 0; i < image.XDim; i++)
                {
                    for (int j = 0; j < image.YDim; j++)
                    {
                        flux.real.data[i, j] = grad.real.data[i, j] * Math.Exp(-(1 / K) * Math.Abs(grad.real.data[i, j]));
                        flux.imag.data[i, j] = grad.imag.data[i, j] * Math.Exp(-(1 / K) * Math.Abs(grad.imag.data[i, j]));
                    }
                }


                // Update image
                for (int i = 0; i < image.XDim - 1; i++)
                {
                    for (int j = 0; j < image.YDim - 1; j++)
                    {
                        image.data[i, j] += lambda * (flux.real.data[i, j] - flux.real.data[i + 1, j]
                            + flux.imag.data[i, j] - flux.imag.data[i, j + 1]);
                    }
                }


                for (int j = 0; j < image.YDim - 1; j++)
                {
                    image.data[image.XDim - 1, j] += lambda * (flux.imag.data[image.XDim - 1, j] - flux.imag.data[image.XDim - 1, j + 1]);
                }

                for (int i = 0; i < image.XDim - 1; i++)
                {
                    image.data[i, image.YDim - 1] += lambda * (flux.imag.data[i, image.YDim - 1] - flux.imag.data[i + 1, image.YDim - 1]);
                }

            }
            public static IntMatrix SegmentAnisotropicGray(System.Drawing.Bitmap Bitmap, int K, int totalIterations, ref System.Drawing.Bitmap GrayDiffused, ref System.Drawing.Bitmap GraySegmented, ref System.Drawing.Bitmap GraySegmented2, ref System.Drawing.Bitmap GrayCurves)
            {
                // Anisotropic diffusion for a gray-scale (single channel) Image.

                // Color image will be automatically converted to brightness...
                RGBIntMatrix imGrayOriginal = RGBIntMatrix.ReadImage(Bitmap);
                DoubleMatrix imGray = imGrayOriginal.Intensity();
                // Update the image in place
                Anisotropic.AnisoDiff(imGray, K, totalIterations);

                GrayDiffused = imGray.WriteBitmapScaled();


                ComplexDoubleMatrix gdGrad = Tools.Gradient(imGray);
                DoubleMatrix f = DoubleMatrix.VectorialAdd(gdGrad);
                DoubleMatrix ctemp = f.Clone();
                DoubleMatrix f3 = Tools.NonMaximaSuppress(ctemp, gdGrad);

                DoubleMatrix t = f.Clone();
                f.Threshold(K);

                IntMatrix labels = Tools.GetEgdes(t, (IntMatrix)f, (IntMatrix)f3);
                RGBIntMatrix save = imGrayOriginal.Clone();
                save.R.MultiplyValue(200);
                save.R.DivideValue(255);
                save.G.MultiplyValue(200);
                save.G.DivideValue(255);
                save.B.MultiplyValue(200);
                save.B.DivideValue(255);

                DoubleMatrix saveInten = save.Intensity();
                IntMatrix curves = new IntMatrix(imGray.XDim, imGray.YDim, 255);

                for (int x = 0; x < imGray.XDim; x++)
                {
                    for (int y = 0; y < imGray.YDim; y++)
                    {
                        if (labels.data[x, y] == 1)
                        {
                            save.r.data[x, y] = 255;
                            save.g.data[x, y] = 255;
                            save.b.data[x, y] = 255;
                            saveInten.data[x, y] = 255;
                            curves.data[x, y] = 0;
                        }
                    }
                }
                GraySegmented = save.WriteBitmap();
                GraySegmented2 = saveInten.WriteBitmapScaled();
                GrayCurves = curves.WriteBitmapScaled();
                return curves;
            }
           
#endregion

# region "For Color Bitmap 2D"
            public static void AnisoDiff(MultiValuedDoubleMatrix image, double K, int iterations)
            {
                for (int i = 0; i < iterations; i++)
                {
                    AnisoDiff(image, K);
                }

            }
            public static void AnisoDiff(MultiValuedDoubleMatrix image, double K)
            {
                ComplexDoubleMatrix[] tmps = new ComplexDoubleMatrix[image.vectorLength];
                for (int i = 0; i < image.vectorLength; i++)
                {
                    tmps[i] = new ComplexDoubleMatrix(image.XDim, image.YDim, 0d);
                }

                ComplexDoubleMatrix grad = new ComplexDoubleMatrix(image.XDim, image.YDim, 0d);

                ComplexDoubleMatrix[] flux = new ComplexDoubleMatrix[image.vectorLength];
                for (int i = 0; i < image.vectorLength; i++)
                {
                    flux[i] = new ComplexDoubleMatrix(image.XDim, image.YDim, 0d);
                }

                // Calculate gradient
                for (int i = 1; i < image.XDim; i++)
                {
                    for (int j = 0; j < image.YDim; j++)
                    {
                        for (int z = 0; z < image.vectorLength; z++)
                        {
                            tmps[z].real.data[i, j] = image.comp[z].data[i - 1, j] - image.comp[z].data[i, j];
                        }
                    }
                }

                for (int i = 0; i < image.XDim; i++)
                {
                    for (int j = 1; j < image.YDim; j++)
                    {
                        for (int z = 0; z < image.vectorLength; z++)
                        {
                            tmps[z].imag.data[i, j] = image.comp[z].data[i, j - 1] - image.comp[z].data[i, j];
                        }
                    }
                }


                DoubleMatrix[] tensor = new DoubleMatrix[3];
                tensor[0] = new DoubleMatrix(image.XDim, image.YDim, 0.0);
                tensor[1] = new DoubleMatrix(image.XDim, image.YDim, 0.0);
                tensor[2] = new DoubleMatrix(image.XDim, image.YDim, 0.0);

                for (int i = 0; i < image.XDim; i++)
                {
                    for (int j = 0; j < image.YDim; j++)
                    {
                        for (int z = 0; z < image.VectorLength; z++)
                        {
                            tensor[0].data[i, j] += tmps[z].Real.data[i, j] * tmps[z].Imag.data[i, j];
                            tensor[1].data[i, j] += tmps[z].Real.data[i, j] * tmps[z].Real.data[i, j];
                            tensor[2].data[i, j] += tmps[z].Imag.data[i, j] * tmps[z].Imag.data[i, j];
                        }


                        double a = tensor[1].data[i, j];
                        double b = tensor[0].data[i, j];
                        double c = tensor[2].data[i, j];

                        double eig = ((a + c) + Math.Sqrt((a - c) * (a - c) + 4 * b * b)) / 2;
                        double theta = Math.Atan2(eig - a, b);
                        grad.real.data[i, j] = Math.Sqrt(eig) * Math.Cos(theta);
                        grad.imag.data[i, j] = Math.Sqrt(eig) * Math.Sin(theta);
                    }
                }


                // Calculate flux
                for (int i = 0; i < image.XDim; i++)
                {
                    for (int j = 0; j < image.YDim; j++)
                    {
                        for (int z = 0; z < image.VectorLength; z++)
                        {
                            flux[z].real.data[i, j] = tmps[z].real.data[i, j] * Math.Exp(-(1 / K) * Math.Abs(grad.real.data[i, j]));
                            flux[z].imag.data[i, j] = tmps[z].imag.data[i, j] * Math.Exp(-(1 / K) * Math.Abs(grad.imag.data[i, j]));
                        }
                    }
                }


                // Update image
                for (int i = 0; i < image.XDim - 1; i++)
                {
                    for (int j = 0; j < image.YDim - 1; j++)
                    {
                        for (int z = 0; z < image.VectorLength; z++)
                        {
                            image.comp[z].data[i, j] += lambda * (flux[z].real.data[i, j] - flux[z].real.data[i + 1, j]
                                + flux[z].imag.data[i, j] - flux[z].imag.data[i, j + 1]);
                        }
                    }
                }


                for (int j = 0; j < image.YDim - 1; j++)
                {
                    for (int z = 0; z < image.VectorLength; z++)
                    {
                        image.comp[z].data[image.XDim - 1, j] += lambda * (flux[z].imag.data[image.XDim - 1, j] - flux[z].imag.data[image.XDim - 1, j + 1]);
                    }
                }

                for (int i = 0; i < image.XDim - 1; i++)
                {
                    for (int z = 0; z < image.VectorLength; z++)
                    {
                        image.comp[z].data[i, image.YDim - 1] += lambda * (flux[z].imag.data[i, image.YDim - 1] - flux[z].imag.data[i + 1, image.YDim - 1]);
                    }
                }

            }


            public static IntMatrix SegmentAnisotropicColor(System.Drawing.Bitmap imname, int K, int totalIterations)
            {
                // Anisotropic diffusion for a color 
                // (multi channel - can be of more than 3 dimensions) Image.

                RGBIntMatrix imOriginal = RGBIntMatrix.ReadImage(imname);
                RGBDoubleMatrix im = (RGBDoubleMatrix)imOriginal.Clone();
                // Convert to Lab color space.
                im = ColorSpaces.RGB2Lab(im);
                // Update the image in place (this time in color space)
                Anisotropic.AnisoDiff((MultiValuedDoubleMatrix)im, K, totalIterations);

                ColorSpaces.Lab2RGB(im).WriteImageScaled(String.Format("ColorDiffused.{0}.png", K));

                im.ImageScale();
                ComplexDoubleMatrix gdGradColor = Tools.GradientMulti(im);
                DoubleMatrix fColor = DoubleMatrix.VectorialAdd(gdGradColor);
                DoubleMatrix ctempColor = fColor.Clone();
                DoubleMatrix f3Color = Tools.NonMaximaSuppress(ctempColor, gdGradColor);

                DoubleMatrix tColor = fColor.Clone();
                fColor.Threshold(K / 3d);

                IntMatrix labelsColor = Tools.GetEgdes(tColor, (IntMatrix)fColor, (IntMatrix)f3Color);
                RGBIntMatrix saveColor = imOriginal.Clone();
                saveColor.R.MultiplyValue(200);
                saveColor.R.DivideValue(255);
                saveColor.G.MultiplyValue(200);
                saveColor.G.DivideValue(255);
                saveColor.B.MultiplyValue(200);
                saveColor.B.DivideValue(255);

                IntMatrix curvesColor = new IntMatrix(im.XDim, im.YDim, 255);

                for (int x = 0; x < im.XDim; x++)
                {
                    for (int y = 0; y < im.YDim; y++)
                    {
                        if (labelsColor.data[x, y] == 1)
                        {
                            saveColor.r.data[x, y] = 255;
                            saveColor.g.data[x, y] = 255;
                            saveColor.b.data[x, y] = 255;
                            curvesColor.data[x, y] = 0;
                        }
                    }
                }
                saveColor.WriteImage(String.Format("ColorSegmented.{0}.png", K));
                curvesColor.WriteImage(String.Format("ColorCurves.{0}.png", K));

                return curvesColor;
            }

            
            #endregion

# region "MetaView 2D"
            public static void AnisotropicDiffuse(ref Matrix3DInt16 Data, int K, int Iterations)
            {
                Int16 Max, Min;
                Max = Data.Maximum();
                Min = Data.Minimum();
                int colorvalue;
                for (int z = 0; z < Data.z; z++)
                {
                    //for each transverse slice
                    System.Drawing.Bitmap B = new System.Drawing.Bitmap(Data.x, Data.y);
                    System.Drawing.Bitmap Bnew = new System.Drawing.Bitmap(Data.x, Data.y);
                    System.Drawing.Bitmap Bseg = new System.Drawing.Bitmap(Data.x, Data.y);
                    System.Drawing.Bitmap Bseg2 = new System.Drawing.Bitmap(Data.x, Data.y);
                    System.Drawing.Bitmap Bcur = new System.Drawing.Bitmap(Data.x, Data.y);
                    for (int x = 0; x < Data.x; x++)
                    {
                        for (int y = 0; y < Data.y; y++)
                        {
                            colorvalue = (int)((Data[x, y, z] - Min) * 255 / (Max - Min));
                            B.SetPixel(x, y, System.Drawing.Color.FromArgb(colorvalue, colorvalue, colorvalue));
                        }
                    }
                    IntMatrix temp = SegmentAnisotropicGray(B, K, Iterations, ref Bnew, ref Bseg, ref Bseg2, ref Bcur);
                    for (int x = 0; x < Data.x; x++)
                    {
                        for (int y = 0; y < Data.y; y++)
                        {
                            colorvalue = Bnew.GetPixel(x, y).R;
                            Data[x, y, z] = (short)((colorvalue * (Max - Min) / 255) + Min);
                            B.SetPixel(x, y, System.Drawing.Color.FromArgb(colorvalue, colorvalue, colorvalue));
                        }
                    }


                }

            }
            #endregion

            # region "MetaView 3D"

            public void AnisoDiff3D(Matrix3DInt16 image, double K,int iterations)
            {
                Matrix3DSingle temp = new Matrix3DSingle(image.x, image.y, image.z);
                float max, min;
                max = image.Maximum();
                min = image.Minimum();
                //make a temporary copy
                for (int i = 0; i < image.x; i++)
                    for (int j = 0; j < image.y; j++)
                        for (int k = 0; k < image.z; k++)
                            temp[i, j, k] = (Single)((image[i, j, k] - min) * 255 / (max - min));

                for (int i = 0; i < iterations; i++)
                {
                    temp = AnisoDiff3D(temp, K);
                }

                for (int i = 0; i < image.x; i++)
                    for (int j = 0; j < image.y; j++)
                        for (int k = 0; k < image.z; k++)
                            image[i, j, k] = (short)((temp[i, j, k] * (max - min) / 255) + min);

            }
            public void AnisoDiff3D(Matrix3DSingle image, double K, int iterations)
            {
                Matrix3DSingle temp = new Matrix3DSingle(image.x, image.y, image.z);
                float max, min;
                max = image.Maximum();
                min = image.Minimum();
                //make a temporary copy
                for (int i = 0; i < image.x; i++)
                    for (int j = 0; j < image.y; j++)
                        for (int k = 0; k < image.z; k++)
                            temp[i, j, k] = (Single)((image[i, j, k] - min) * 255 / (max - min));

                for (int i = 0; i < iterations; i++)
                {
                    temp = AnisoDiff3D(temp, K);
                }

                for (int i = 0; i < image.x; i++)
                    for (int j = 0; j < image.y; j++)
                        for (int k = 0; k < image.z; k++)
                            image[i, j, k] = (Single)((temp[i, j, k] * (max - min) / 255) + min);

            }

            public Matrix3DSingle AnisoDiff3D(Matrix3DSingle image, double K)
            {
                Matrix4Ddouble Gradients = new Matrix4Ddouble(image.x, image.y, image.z, 6);
                float change;
                // Calculate gradient
                for (int i = 1; i < image.x - 1; i++)
                {
                    for (int j = 1; j < image.y - 1; j++)
                    {
                        for (int k = 1; k < image.z - 1; k++)
                        {
                            Gradients.data[i, j, k, 0] = image[i - 1, j, k] - image[i, j, k];
                            Gradients.data[i, j, k, 1] = image[i + 1, j, k] - image[i, j, k];
                            Gradients.data[i, j, k, 2] = image[i, j + 1, k] - image[i, j, k];
                            Gradients.data[i, j, k, 3] = image[i, j - 1, k] - image[i, j, k];
                            Gradients.data[i, j, k, 4] = image[i, j, k + 1] - image[i, j, k];
                            Gradients.data[i, j, k, 5] = image[i, j, k - 1] - image[i, j, k];
                            //         global::System.Windows.Forms.MessageBox.Show(image[i,j,k].ToString());
                            change = (float)(0.16 * (
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 0] / K) * (Gradients.data[i, j, k, 0] / K))) * Gradients.data[i, j, k, 0] +
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 1] / K) * (Gradients.data[i, j, k, 1] / K))) * Gradients.data[i, j, k, 1] +
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 2] / K) * (Gradients.data[i, j, k, 2] / K))) * Gradients.data[i, j, k, 2] +
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 3] / K) * (Gradients.data[i, j, k, 3] / K))) * Gradients.data[i, j, k, 3] +
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 4] / K) * (Gradients.data[i, j, k, 4] / K))) * Gradients.data[i, j, k, 4] +
                                Math.Exp(-1 * Math.Abs((Gradients.data[i, j, k, 5] / K) * (Gradients.data[i, j, k, 5] / K))) * Gradients.data[i, j, k, 5]));
                            image[i, j, k] += change;

                        }
                    }
                }

                return image;
            }

            #endregion
            
            
        }
    }
}